Polynomial Approximation with Derivatives

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from math import factorial

A Brief Intro to sympy

Here we import sympy, a package for symbolic computation with Python.

import sympy as sp
sp.init_printing()

Next, we make a (symbolic) variable $x$ from which we can then build more complicated expressions:

sp.var("x")
x

Build up an expression with $x$. Assign it to expr. Observe that this expression isn't evaluated--the result of this computation is some Python data that represents the computation:

g = sp.sin(sp.sqrt(x)+2)**2
g

Next, take a derivative, using .diff(x).

g.diff(x)

g.diff(x, 4)

Use .subs(x, ...) and .evalf() to evaluate your expression for $x=1$.

g.subs(x,1)

g.subs(x, 1).evalf()

Helper function that takes symbolic functions as argument and plot them

def plot_sympy(my_f, my_pts, **kwargs):
    f_values = np.array([my_f.subs(x, pt) for pt in my_pts])
    plt.plot(pts, f_values, **kwargs)

Polynomial Approximation

f = 1/(20*x-10)
f

f.diff(x)

f.diff(x,2)

Write out the degree 2 Taylor polynomial about 0:

taylor2 = (
    f.subs(x, 0)
    + f.diff(x).subs(x, 0) * x
    + f.diff(x, 2).subs(x, 0)/2 * x**2
)
taylor2

Plot the exact function f and the taylor approximation taylor2

pts = np.linspace(-0.4,0.4)

plot_sympy(taylor2, pts, label="taylor2")
plot_sympy(f, pts, label="f")
plt.legend(loc="best")
plt.axis('equal')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('function values')

Behavior of the Error

Let's write the taylor approximation for any degree n, and define the error as f - tn

n = 2

tn = 0
for i in range(n+1):
    tn += f.diff(x, i).subs(x, 0)/factorial(i) * x**i

plot_sympy(tn, pts, label="taylor")
plot_sympy(f, pts, label="f")
plt.legend(loc="best")
plt.axis('equal')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('function values')

error = f - tn

plot_sympy(error, pts, label="error")
plt.legend(loc="best")
plt.ylim([-1.3, 1.3])
plt.axis('equal')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('error')

To get a better idea of what happens close to the center, use a log-log plot:

# plot only points close to zero [10^(-3),10^(0.4)]
plt.figure(figsize=(10,6))
pos_pts = 10**np.linspace(-3, 0.4) 
err_values = [abs(error.subs(x, pt)) for pt in pos_pts]
plt.loglog(pos_pts, err_values)
plt.grid()
plt.xlabel("$x$")
plt.ylabel("Error")